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Abstract 

We present results for Charmonium spectroscopy using Non-Relativistic 
QCD (NRQCD). For the NRQCD action the leading order spin-dependent 
and next to leading order spin-independent interactions have been included 
with tadpole-improved coefficients. We use multi-exponential fits to multi- 
ple correlation functions to extract ground and excited S states. Splittings 
between the lowest S, P and D states are given and we have accurate val- 
ues for the S state hyperfine splitting and the Xc fine structure. Agreement 
with experiment is good - the remaining systematic errors are discussed. 

PACS numbers: 12.38.Gc, 14.40.Gx, 14.65.Dw, 12.39.Hg 



1 Introduction 

The study of heavy-heavy mesons is important for Lattice Gauge Theory not only 
because of the availability of experimental data for comparison but also because 
such systems allow a quantitative study of systematic errors which arise in lattice 
simulations at present. To study heavy- heavy mesons we use Non-Relativistic 
QCD (NRQCD) [[l], and previously we have reported a very successful study 
of the Bottomonium system|3|]. This allowed the extraction of two fundamental 
parameters in QCD[|J, the b-quark mass || and the strong coupling constant a s 
[0. Here we report on a similar study of the Charmonium spectrum. 
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The starting point of NRQCD is to expand the original QCD lagrangian in 
powers of v 2 , the typical quark velocity in a bound state. For the J/ty system 
v 2 ~ 0.3. Thus we systematically include relativistic errors order by order in v 2 
away from a No n- Relativistic limit. Our action is the same one as used in Q where 
relativistic corrections OiMv ) have been included. This means that systematic 
errors from relativistic corrections will be 0(M c v 6 ) (= ~ 30 — 40 MeV) for the 
J/ty system i.e. 10% in spin- independent splittings and 30% in spin-dependent 
splittings. This is considerably less accurate than for the T case0 because v 2 
is about a factor of 3 larger here. Other sources of systematic error include 
discretisation errors and errors from the absence of virtual quark loops because 
we use quenched configurations generated with the standard plaquette action. 
Finite volume errors should be negligible because of the relatively small size of 
the J/\I/ system. 

Shown in Figures ([I]) and (|2|) is the spectrum for Charmonium using Lat- 
tice NRQCD. The spectrum was calculated using an ensemble of 273 gauge field 
configurations generated with the standard Wilson action at j3 — 5.7J7J. To set 
the scale we fix our simulation result for the spin-averaged IP-IS splitting to its 
experimental value of 458 MeV. This gives a~ l = 1.23(4) GeV, where the uncer- 
tainty is purely statistical. Since we are working in the quenched approximation 
this value can be and is different both from that obtained at the same value of (3 
using light hadron spectroscopy || or using Upsilon spectroscopy^]. We expect 
a value fixed from heavyonium to be more accurate than that from light hadron 
spectroscopy because spin-independent splittings in the heavy quark sector are 
independent of quark mass to a good approximation and systematic errors are 
under better control |T]|. 

To fix the bare quark mass in the action, M c °, we plot a dispersion relation 
correct up to 0(v A ) for the r] c . M c ° is then tuned until the simulation value for 
the kinetic mass is equal to the experimental value of the mass of the r] c (2.98 
GeV). We find that using aM c °=0.8 gives M(r/ C )=3.0(l) GeV with a" 1 = 1.23(4) 
GeV. 

In Figure ([I]) the whole Charmonium spectrum is shown and in Figure (|2|) 
the spin-dependent splittings are shown in more detail. In Figure (^) it can 
be seen that although the general pattern of splittings for the S and P states 
is reproduced well, systematic errors are visible above the statistical errors. It 
should then be possible in the future to observe systematic improvements to the 
current calculation, when higher order relativistic corrections are included and 
further discretisation and quenching errors are removed. 

We give details in section 2 of our evolution equation and the quark Greens 
function used to make up meson correlation functions. Section 3 describes the 
results from the simulation using multi-exponential fits. We illustrate the need 
for multiple smearing functions to obtain smaller statistical errors. Section 4 
compares simulation results to experiment and section 5 contains our conclusion. 
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Figure 1: NRQCD simulation results for the spectrum of the J/*f> system using an 
inverse lattice spacing of 1.23 GeV, fixed from the spin-averaged IP-IS splitting. 
The 1 S'o mass is fixed at 3.0 GeV, from a fit to the kinetic mass. Experimental 
values are indicated by dashed lines. Error bars are shown where visible, and 
only indicate statistical uncertainties. 

2 Evolution Equation and Quark propagators 

One of the advantages of the formulation of NRQCD is that it involves a simple 
difference equation in the temporal direction. This allows the evolution of the 
quark Green function as an initial value problem which can be solved with one 
sweep through the lattice. We define our quark Green function to be initially 
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and then continue to evolve using 
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On the lattice, the kinetic energy operator is 
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Figure 2: Simulation results for the spin structure of the J/^f family, using an 
inverse lattice spacing of 1.23 GeV. The energies of the spin- averaged S and P 
states have been set to zero. Error bars for points are statistical. 
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The first two terms in SH are spin-independent relativistic corrections and the 
next two are spin- dependent correction terms which contribute to the P and S 
spin splittings respectively. The last two terms come from finite lattice spacing 
corrections to the lattice Laplacian and the lattice time derivative. A is the 
symmetric lattice derivative, A^ 2 ^ is the lattice form of the Laplacian and A^ 4 ) is 
a lattice version of the continuum operator Df. We used the standard traceless 
cloverleaf operators for the chromo-electric and magnetic fields, E and B. The 
parameter n is introduced to remove instabilities in the heavy quark propagator 
caused by the highest momentum modes of the theory Q. For our simulations at 
(3 = 5.7 and with a bare mass for the c quark in lattice units of 0.8, we set n = 4. 

The coupling constants appearing in equation (|J) can be calculated by 
matching NRQCD to full QCD || |10j. At tree level all the coefficients are one. 



The largest radiative corrections are believed to be tadpole contributions pl| . We 
take care of these by using the method suggested in [TT| where all the U's are 
redefined by 

U,(x) - ^ (5) 
with m the fourth root of the plaquette (at (3=5.7 we use uq = 0.861). Since the 
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cloverleaf expression involves the evaluation of a plaquette this renormalization 
will have the effect of redefining E and B via 



(«> 

which will strongly affect spin-dependent splittings. With the dominant tadpole 
contributions thus removed, we use the tree level values for the q's. The only 
remaining free parameters are the bare quark mass M c ° and the bare coupling 
constant g which appear in the original QCD Lagrangian. All the details of the 
quark evolution up to this point are identical to those in ||. In the following 
some of the technical details differ slightly. 

Given the quark propagators in equation (|2|) it is relatively straightforward to 
combine them appropriately to form meson propagators with specific quantum 
numbers. Using the notation of || we take ip> to create a heavy quark and x 
to create a heavy anti-quark. Then the following interpolating operator creates 
a meson centred on the point X\ : 

^^ f (fi)r(fi -X2)x\x 2 ). (7) 

X2 

Local meson operators are tabulated in ||. Here we generalise the operators 
to include 'smearing functions'. For S states the meson operator F becomes 
Q 4>{xx —X2) where Q is a 2 x 2 matrix in spin space giving the quantum numbers 
of the meson and is a simple approximation to the wavefunction. For P states, 
</> also becomes a p wavefunction, which can be thought of as the derivative of 
a spherically symmetric function ||. In general T is a sum of spin matrices 
multiplying different smearing functions, generalising the operators in ||. For 
the wavefunctions <fi we use here wavefunctions from a 1/r potential with their 
spread adjusted to match the size of the appropriate meson. 
For meson propagators at zero momentum we then have 



yi.3/2 



G meson (p = 0, t) = £ Tr G\(y 2 ) r^j/i - V2) G t (yi) 



(8) 



with 



G t (y) = -£G t (y-x)T( sc \x). (9) 



r sc (a;) and T sk (x) refers to the meson operator T(x) = Q<p(x) with the smearing 
function <p{x) at the source or sink respectively and enumerated by the integer n sc 
or risk, n — 1 corresponds to the ground state meson, n = 2 to the first radially 
excited state. G t is obtained using equations (1) and (2) with Sgp — ► T^ sc \x). The 
trace is over color and spin. The convolutions are evaluated using Fast Fourier 
Transforms. 
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We also study finite momentum propagators for the 1 So meson, given by: 

G meson {p,t) =J2 Tr [G\{yi)nG t {y x )\ e-*** (10) 



yi 



Using the notation 2S+1 Lj, we have looked at meson propagators for the 
following states: 1 S , 3 Si, 1 P i , 3 P , 3 Pi, 3 P 2 for both the E and T representation 
and the l D 2 in the T representation. For the S states, smearing functions both for 
the ground and first radially excited state were used as well as a local 5 function 
(n = loc). From this all possible combinations of smearing at the source and sink 
were formed making a 3 x 3 matrix of correlation functions. For the P and D states 
only the ground state smearing function was used at the source. We calculated 
the dispersion relation for the 1 Sq by looking at the meson propagator for small 
momentum components using (n sc ,n s k) = (loc, loc) and (l,loc). To maximize 
our statistics we use all color and spin indices at the source when calculating 
our meson propagators. For the 3 Si, l P\, 3 Pi, 3 P 2 and 1 D 2 we average over 
polarization directions making a total of 30 S, P and D meson propagators to 
analyze. 



3 Simulation results 

In the simulation we used 273 quenched gluon field configurations on a 12 3 x24 
lattice at /3 = 5.7 generously supplied by the UKQCD collaboration |7| . They were 
fixed to Coulomb gauge using a Fourier accelerated steepest descents algorithm 
[p"2]j with a cutoff on [d ■ A} 2 of 10~ 6 . Due to the relatively small size of the 
J/ty it is possible to use more than one starting site on a spatial slice. We also 
use more than one starting point in time to increase statistics. In this case we 
used 8 different spatial origins and 2 different starting times at timeslice 1 and 
12. If we bin the spatial origins together we find significant correlation, whereas 
binning together two propagators with an initial timeslice of 1 and 12 but with 
the same spatial origin gives little or no correlation at all. For most of our fits 
we bin together all the correlation functions from a given configuration, except 
when doing multiple-exponential multiple- correlation fits for the 1 S'o and 3 Si case. 
Here we only bin on spatial origin and having the increased sample size from the 
time direction significantly improves the fit. We also checked, however, that 
fitting with all data unbinned produces a worse \ 2 than when all data is binned, 
another indicator of spatial correlations. 

In NRQCD, as in QCD, there are two free parameters, the bare coupling 
constant g and the bare quark mass M c °. We fix g implicitly when we set the 
scale a -1 . To fix M° we tune so that the simulation result for the kinetic mass 
of the 1 S'o agrees with the experimental value of the mass of the rj c (2.98 GeV). 



6 



For this we find E-p for several different momenta of the 1 5'o and fit to the form 
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simultaneously for momenta components (1,0,0), (1,1,0), (1,1,1), (2,0,0) in units 
of 2ir/12a. M k i n is taken to be the rest mass of the r] c . C\ should take the value 
1.0 in a fully Lorentz invariant theory. Instead we find the value 1.7(1) - this is 
because of relativistic corrections that have not been included. The mass in the 
P 4 term differs from that in the P 2 term by the cube root of 1.7 i.e. 20%. Since 
we have included no relativistic corrections to the P 4 term we would expect it 
to be correct only to leading order i.e. 30%, and the difference we observe is 
consistent with that. We expect results closer to 1.0 for the T because this is 
a more non-relativistic system. Indeed there || we find a tendency for C\ to 
be larger than 1.0 but consistent with 1.0 within rather larger errors of size 0.3. 
The last term in equation (|ll]) is a non-rotationally-invariant term allowed on 
the lattice but C<i is found to be —0.1(1) consistent with zero. This indicates 
that no discretisation errors are visible in the dispersion relation once the 0(a 2 ) 
terms in the heavy quark action have been taken care of in equation (^). A fit 
with the extra P 6 relativistic correction in was tried and no significant signal for 
it was found. Conversely a fit with just the leading order P 2 term in was tried 
but gave a very poor Q value. This suggests that with the particular momentum 
components used a fit including terms up to P 4 is appropriate. Using a bare 
quark mass of oM c ° = 0.8 gives M km a = 2.429(7) or M kin = 3.0(1) GeV using 
a -1 = 1.23(4) GeV. All simulation results quoted here are from using this value of 
the bare quark mass. The error on the bare quark mass is then of order 10% from 
both statistical errors in a -1 and systematic errors from higher order relativistic 
corrections. 

In Figures @ and (U) we show effective masses for the 1 So and 1 P\ states 
respectively. We use the naive definition m e ff(t) = —log(G(t+ 1)/G{t)) together 
with bootstrap errors. From the S state plots it is clear that smearing has the 
effect of producing an earlier plateau in the effective mass. Although the statis- 
tical errors have increased for the smeared cases as compared to the local-local 
case the earlier plateau allows fitting to take place closer to the origin and ul- 
timately produces better errors. For the first excited state a plateau cannot be 
seen for the effective mass and the signal ultimately decays to the ground state. 
A better transient plateau was seen for the excited S state in the T spectrum 
at P = 6.0 This reflects the fact that at higher (3 values the excited states 
have smaller masses in lattice units and last for longer times. For the P state the 
signal/noise ratio is much poorer than that for the S state, as expected. 
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Figure 3: 1 Sq Effective masses by (source, sink). 
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3.1 Fitting Results for the 1 Sq and 3 Si and the singlet P 
and D states 

We use a variety of fitting routines to extract high precision ground state masses 
for the 1 Sq and 3 S\ as well as masses for their first radially excited states. We 
have used in general the same fitting procedures which are described in more 
detail in ||. 

Multi-exponential fits allow a fit to the correlation function at much earlier 
times than single-exponential fits, thus reducing the noise. A fit to n exponentials 
allows confidence in the masses of the first n — 1 states. Since, as described above, 
excited states die very rapidly at low (3, it is much harder to get a value for an 
excited state mass at (3=5.7 than at (3 = 6.0. This is reflected in our errors. It is 
also true, however, that the ground state plateau appears earlier and the use of 
many exponentials to get to early times is not as important at [3 = 5.7 as at (3 = 
6.0. 

The first type of fit we do is that to a matrix of correlation functions: 

Nexp 

G meS on(n sc ,n sk ]t) = flKc^'Kt,^^' 1 (12) 

k=l 

For the S states we use the combination n sc = 1,2 and n s k = 1,2 forming a 2 
x 2 matrix. Then we perform fits for N exp = 2 and 3. Our fitting procedure 
inverts the covariance matrix using the svd algorithm. We have sufficiently good 
statistics that we are able to keep all eigenvectors of the covariance matrix and 
achieve a good fit [|l3]| . 



For the second fit a row of correlation functions is formed and fitted to 

G m eson(n sc ,loC]t) = b(n sc , k) 6~ Ek t (13) 

k=l 

We use the correlation functions (n sc ,n s k) with n sc = 1,2 n s k = loc. Again fits 
use N exp = 2, 3. 

In Tables [I] and [5] are results from the row and matrix fits for the 1 S'o and 3 S\. 
The errors stated are those causing a change 8\ 2 = 1 and we also quote the quality 
of the fit, Q. For an acceptable fit Q should be in the range 0.01 to 0.9 and ideally 
Q > 0.1. To improve our statistics we only bin correlation functions which start 
from different spatial origins but not ones which have different starting timeslices. 
This has little effect on the central value but does increase the Q value giving us 
more confidence in the fit. 

From both tables it is clear that an accurate ground state mass can be obtained 
at very early times. Only a t min of 2 gives an unacceptable Q for the 2 exponential 
fit. Adding a 3rd exponential produces an acceptable fit, although we don't take 
this value because Q increases further as t m i n is increased. This contrasts with 
the higher t m i n needed for T spectroscopy at (3 = 6.00. The masses we obtain 
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Figure 4: x Pi Effective masses by (source, sink). 



are independent of the type of fitting routine within errors, although the values 
for Q are lower for the matrix fits. At this point it is constructive to test how 
effective the multiple exponential fits are for the ground states at /3 — 5.7. In 
Table || are values for a single exponential fit to the (n sc , n s k) = (1, loc) and (1, 1) 
for the 1 S state. In both cases an acceptable Q requires t min of 6, significantly 
larger than for the multiple exponential fit. We choose fitted values 0.6182(7) 
and 0.697(1) for the 1 S and 3 Si ground states respectively. 

For the first excited state the choice of fitted value is far more difficult. To 
have confidence in the value we should use a 3 exponential fit although this gives 
larger errors in the fitted masses. We look for both a steady value in the fitted 
mass as t min is changed and a steady value for Q. It is also useful to look at the 
amplitude for the second excited state in the 3 exponential fit to see at what t min 
values it has decayed away. 

For the 1 So row fit we choose a value 1.17(5) for the excited state mass (average 
of t m i n = 3,4,6) and from the matrix fit 1.18(4) (average of t m i n = 3,4,5). There 
is then agreement within errors between the two fits and we choose 1.17(5) as 
the global average. For the 3 Si state there is a significant deterioration in the 
Q values over those for the 1 Sq and the fitting errors are slightly larger. This is 
presumably a reflection of the additional noise in the 3 Si channel coming from 
the 1 S . For the row fit a value of 1.19(7) (average for t min = 4,5,6) is chosen and 
a value of 1.22(3) (average for t min =3,4,5) from the matrix fit. A global average 
for the excited 3 Si is chosen to be 1.20(7). All the fitted values are collected in 
Table 0. 

In Tables f| and |5| are the amplitudes from the various fits for particular 
values of t m i n /t max . The value of t m i n /t max used was that where the fit for the 
first excited state was closest to the average result quoted above. In both row 
and matrix fits it was found that the amplitude for a second excited state {k = 
3) is essentially zero. This indicates that contamination from higher states in our 
fits is negligible. From the amplitude results we can see that n sc =l has strongest 
overlap with the ground state and n sc = 2 has strongest overlap with the first 
excited state, as planned. Thus our smearing functions are projecting out the 
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Q 


fits to (l,loc) 


2 


2/24 


0.6171(6) 


1.172(6) 




2 x 10~ a 


and (2,loc) 




3/24 


0.6178(6) 


1.16(1) 




0.65 






4/24 


0.6176(6) 


1.16(1) 




0.64 






5/24 


0.6179(7) 


1.14(1) 




0.79 






6/24 


0.6182(7) 


1.21(5) 




0.94 






7/24 


0.6183(7) 


1.27(8) 




0.93 




3 


2/24 


0.6180(7) 


1.15(2) 


1.8(6) 


0.38 






3/24 


0.6177(20) 


1.15(4) 


1.8 ± 15 


0.53 






4/24 


0.6181(6) 


1.16(2) 


1.8(1.2) 


0.79 






5/24 


0.6183(7) 


1.30(16) 


1.7(6) 


0.94 






6/24 


0.6183(7) 


1.19(8) 


1.8(5) 


0.87 






7/24 


0.6183(7) 


1.25(24) 


1.8(8) 


0.85 


fits to 


2 


3/24 


0.6185(6) 


1.18(2) 




0.06 


(1.1), (1,2) 




4/24 


0.6183(6) 


1.17(3) 




0.15 


(2,1), (2,2) 




5/24 


0.6178(6) 


1.16(4) 




0.25 






6/24 


0.6177(6) 


1.08(6) 




0.16 






7/24 


0.6181(6) 


0.90(6) 




0.42 




3 


3/24 


0.6180(6) 


1.19(2) 


1.6(5) 


0.27 






4/24 


0.6178(6) 


1.14(4) 


2.1(6) 


0.23 






5/24 


0.6179(6) 


1.21(7) 


1.7(6) 


0.16 






6/24 


0.6180(6) 


1.26(11) 


2(1) 


0.18 






7/24 


0.6181(6) 


0.91(6) 


2(1) 


0.33 



Table 1: Examples of simultaneous multi-exponential fits to the 1 Sq using row 
and matrix fits respectively. 

required state and suppressing the others, although our smearing functions are 
clearly not optimal. It may be better to use the output wavefunctions to produce 
input smearing functions in an improved calculation. To illustrate the quality of 
the multi-exponential fits into early times we have plotted in Figure 5 effective 
amplitude plots with the fitted parameters quoted in Tables |] and 

For the P and D states multiple exponential fits are not possible because we 
have included only the ground state smearing function in the simulation. Instead 
a single exponential fit was performed to the (n sc , n s k) = (1, 1) meson propagators 
of the x Pi and 1 D 2 . The results are shown in Tables ^| and 0. Reasonable errors 
are obtained at t min values of 6 where single exponential fits were acceptable for 
the S states. Ratio fits were also done to the 1 S in both cases but the results 
and errors remained the same showing there is no correlation between these states 
and the 1 5'o. To isolate the ground state early on and achieve better errors higher 
radial smearing functions need to be added. Work has begun on this for the x Pi 
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aE x 
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fits to (l,loc) 


1 


5/24 


0.6188(8) 


0.01 






6/24 


0.6184(8) 


0.66 
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0.6183(8) 


0.72 


fits to (1,1) 
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4/24 


0.6184(8) 


0.05 






5/24 


0.6181(8) 


0.22 






6/24 


0.6181(8) 


0.18 






7/24 


0.6182(8) 


0.15 



Table 2: Examples of single exponential fits to the 1 S , . 





N 

1 ' exp 


t rain /t max 




aE 2 


aE 3 


Q 


fits to (l,loc) 


2 


2/24 


0.6951(8) 


1.247(7) 




4 x 10~ 5 


and (2,loc) 




3/24 


0.6961(8) 


1.23(1) 




0.23 






4/24 


0.6958(9) 


1.22(2) 




0.23 






5/24 


0.6961(9) 


1.18(2) 




0.46 






6/24 


0.6966(9) 


1.21(5) 




0.56 






7/24 


0.6968(10) 


1.25(8) 




0.56 




3 


2/24 


0.6964(9) 


1.21(4) 


1.9(9) 


0.10 






3/24 


0.6957(9) 


1.20(4) 


1.9(1.4) 


0.17 






4/24 


0.6964(10) 


1.16(5) 


1.9(1.3) 


0.47 






5/24 


0.6967(10) 


1.22(8) 


1.9(5) 


0.55 






6/24 


0.6966(7) 


1.19(6) 


1.9(3) 


0.41 






7/24 


0.6969(10) 


1.25(16) 


1.9(2) 


0.40 


fits to 


2 


3/24 


0.6970(8) 


1.22(1) 




0.04 


(1,1), (1,2) 




4/24 


0.6967(8) 


1.21(3)) 




0.05 


(2,1), (2,2) 




5/24 


0.6965(8) 


1.24(5) 




0.07 






6/24 


0.6966(8) 


1.31(9) 




0.09 






7/24 


0.6967(9) 


0.95(8) 




0.08 




3 


3/24 


0.6966(8) 


1.23(2) 


1.7(6) 


0.08 






4/24 


0.6965(8) 


1.20(3) 


1.8(6) 


0.06 






5/24 


0.6964(8) 


1.23(4) 


2.0(2.8) 


0.04 






6/24 


0.6969(8) 


1.46(13) 


1.8(1.3) 


0.06 






7/24 


0.6967(9) 


1.00(9) 


1.9(1.3) 


0.07 



Table 3: Examples of simultaneous multi-exponential fits to the 3 Si using row 
and matrix fits respectively. 
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Fit 


tmin / tmax 


k 


a(n SCjSfc = 1, k) 


a{n SCySk = 2, k) 


N eX p 2 


4/24 


1 


0.681(1) 


-0.1188(8) 


for % 




2 


0.18(9) 


0.52(2) 


N eX p 2 


5/24 


1 


0.700(3) 


-0.164(1) 


for 3 S l 




2 


0.29(2) 


0.53(5) 



Table 4: Examples of fit results for amplitudes a(n SCjS fc, k) 



Fit 


^min / ^"max 


k 


b(n sc = l,k) 


b(n sc = 2, k) 


N e xp = 2 


4/24 


1 


0.1037(7) 


-0.0184(4) 


for % 




2 


0.032(3) 


0.064(2) 


N eX p 2 


5/24 


1 


0.103(1) 


-0.0253(4) 


for 3 S l 




2 


0.036(7) 


0.069(3) 



Table 5: Examples of fit results for amplitudes b(n sc , k) 



state. 

3.2 Fits to Spin Splittings 

As described earlier, spin splittings are very dependent on the tadpole improved 
coupling constants q. This makes the spin-splittings a good test of the tadpole- 
improvement scheme. It is also true that potential models find it hard to produce 
spin-splittings in agreement with experiment so we would hope that they are also 
a good test of the differences between a full calculation in QCD, such as ours, 
and a potential model. 

Since meson correlation functions of given / from the same configuration are 
highly correlated we produce a bootstrap ensemble of ratios of correlation func- 





N eX p tmin /tmax 


aEi 


Q 


fits to (1,1) 


1 3/24 


1.059(4) 


0.45 




4/24 


1.052(5) 


0.68 




5/24 


1.049(7) 


0.66 




6/24 


1.046(9) 


0.62 




7/24 


1.048(14) 


0.55 



Table 6: Example of a 1 Pi fit. 
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N eX p train /^max 


aE x 


Q 


fits to (1,1) 


1 3/24 


1.35(1) 


0.62 




4/24 


1.32(2) 


0.77 




5/24 


1.30(3) 


0.78 




6/24 


1.26(5) 


0.78 




7/24 


1.26(9) 


0.72 



Table 7: Example of a 1 D 2 fit. 
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Figure 5: 1 S Effective amplitudes G(t) -e^ 1 '* from two-exponential row fits (1,£) 
(2,£) and two-exponential matrix fits (1,1) (2,2) with t min = 4:,t max = 24 . 
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Splitting 






a5E 


Q 


3 Q IP 


1 


4/24 


0.0794(3) 


4.0 x 10" 5 






6/24 


0.0784(4) 


0.35 






8/24 


0.0784(4) 


0.32 






10/24 


0.0783(5) 


0.21 






12/24 


0.0778(6) 


0.25 


3 p 3 p 
r 2E — -M) 


1 


3/13 


0.090(2) 


0.94 






4/13 


0.089(4) 


0.91 






5/13 


0.090(6) 


0.85 






6/13 


0.086(9) 


0.81 




1 


3/13 


0.045(1) 


0.99 






4/13 


0.046(3) 


0.99 






5/13 


0.045(4) 


0.99 






6/13 


0.044(6) 


0.97 



Table 8: Examples of ratio fits for spin-splittings 



tions to find spin splittings. From this we fit to a single exponential 

Ratioit) = Ae~ SEt (14) 

We use correlation functions with (n sc ,n s k) = (1,1) and bin on time and spatial 
origin. We find very high Q values in general. Shown in Table § are values 
obtained for various combinations of spin-splittings using equation (|T4|). The SE 
obtained for the 3 Si to 1 S'o ratio fit is in agreement with that obtained from the 
separate row and matrix fits of Tables |TJ and |3|. To estimate 5E for higher radial 
excitations we have used a correlated 5E fit. This is a fit to the form 



G m eson Aij^sc) loC\ t) — Ca(^sc ; 6 

k=l 

*y, exp 

G mesonB (n sc Joc;t) = c B (n sc ,l) e~^ +SE > t + ]T c B (n sc , k) e~ E * 4 (15) 

k=2 

with n sc = 1,2 for each meson. The results shown in Table |9| show that the 
3 S\ — 1 Sq splitting can be obtained at early times with smaller errors than in the 
ratio fit. Presumably extra excited states have been absorbed in the extra terms 
in the correlated fit. We are unable to obtain a clear signal for a 2S hyperfine 
splitting although the correlated 5E fit above and the individual matrix fits give 
an indication of such a splitting at early times. 
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train/ tmax 


1% 


2% 


l 3 5i -l 1 ^ 




Q 


3/24 


0.6179(6) 


1.17(1) 


0.0779(3) 


1.23(1) 


0.29 


4/24 


0.6178(6) 


1.17(1) 


0.0778(3) 


1.24(2) 


0.33 


5/24 


0.6180(6) 


1.16(2) 


0.0777(3) 


1.19(2) 


0.72 


6/24 


0.6183(6) 


1.20(4) 


0.0780(4) 


1.20(4) 


0.88 


7/24 


0.6183(7) 


1.20(6) 


0.0781(4) 


1.20(6) 


0.82 


8/24 


0.6184(7) 


1.16(11) 


0.0781(4) 


1.25(13) 


0.76 



Table 9: Example of correlated SE fit for the 3 Si and 1 S states 





Simulation Results 




0.6182(7) 




0.697(1) 


2^0 


1.17(5) 


2 3 5i 


1.20(7) 


l 1 ^ 


1.05(1) 




1.30(4) 


d Si — 1 So 


0.0782(4) 


3 P 2 - 3 Po 


0.088(8) 


3 P 2 - 3 Pi 


0.044(5) 


3 Pi - 3 Po 


0.044(3) 


3 PcM — 1 Pi 


0.010(1) 



Table 10: Fitted dimensionless energies. 
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4 Comparison with Experiment 



In Table |10| we give the dimensionless splittings obtained from our fitting proce- 
dure. To compare simulation results to experiment it is necessary to fix the scale 
a -1 . We choose the spin-averaged IP-IS splitting to do this. By spin-averaged 
splitting we mean the splitting between spin-averaged states. The spin-averaged 
S state has mass 0.25 x [3m( 3 S'i) + m( 1 So)]. The spin-averaged P state has either 
the mass of the 1 P\ or mass 1/9 x [5m( 3 P2) + 3m( 3 Pi) + m( 3 Po)]. These two P 
masses are the same in potential models and experimentally they do seem to be 
very close although the mass of the 1 P\ needs confirmation ||14|| . In our simulation 
the two masses are slightly different (see Table [10]). We will use m( 1 Pi) because 
of the previously noted disagreement with experiment in the P fine structure. 
The difference in value of the a _1 s obtained is within the statistical error. 

The spin-averaged IP-IS splitting has the advantage of being independent 
of any errors in spin-dependent terms and of being experimentally known to be 
independent of the heavy quark mass in the c, b, region. This gives much less sys- 
tematic uncertainty than, for example, in light hadron spectrum determinations 
of a -1 . In the T spectrum calculation J| it was possible to see a difference in a' 1 
between that fixed from the 2S-1S splitting and that fixed from the IP-IS split- 
ting. Here both our statistical error on the 2S state and our expected systematic 
error from relativistic corrections are too large for this to be possible. 

Using the values in Table |TUj we find a" 1 = 1.23(4) GeV from the 1P( 1 Pi)-1S 
splitting. In Table |TT| we compare the splittings obtained from this simulation 
with experimental results. The results are plotted in Figures ([!]) and ©. It is 
important to remember that there is a potential 30 — 40 MeV systematic error in 
all splittings coming from relativistic corrections not included in the heavy quark 
action. Table [11] and the figures do not include the statistical error in a -1 in 
their quoted errors since all the splittings are correlated. Table [TT] does, however, 
include this error for the hyperfine splitting since this is very sensitive to shifts 
in the bare quark mass allowed by uncertainties in a -1 (the hyperfine splitting 
behaves as 1/Mq in perturbation theory, see equation (16) below). Using the \c 
average for IP would give a~ l = 1.20, at the lower end of the range for a -1 from 
the 1 P 1 . 

As discussed earlier, the statistical error on the 2S state is too large to see any 
significance in the fact that it is slightly higher than experiment. The direction of 
the slight disagreement is the same as that for the T spectrum || . There it seems 
clear that the correction of 0(a 2 ) errors in the gluon action and unquenching will 
produce agreement with experiment || |1(J. To test this for the ^ we will need 



to reduce the statistical errors and systematic errors from the heavy quark action 
in the 2S state. 

The expected shift in the IS state from gluonic 0(a 2 ) effects is 0.006 in lattice 
units for this simulation. This is calculated either perturbatively from the wave- 
function at the origin || or non-perturbatively using a lattice potential model 
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Simulation Results [GeV] 


Experiment [GeVl 

1 L J 


2 L S - 


1% 


0.68(6) 




2 3 S 1 - 




0.62(8) 


0.589(1) 


l D 2 - 


5o 


0.84(5) 




3 D, - 


So 




0.791(3) 


*S! ~ 


% 


0.096(2) 


0.118(2) 


3 P 2 - 


3 A> 


0.11(1) 


0.141(1) 


3 P 2 " 




0.054(6) 


0.0456(1) 


3 PcM ~ 




0.012(1) 


0.0008(3)* 



Table 11: NRQCD spectrum results and comparison with experiment for a 
1.23 GeV and aM c ° = 0.8. * requires confirmation. 



15(1 - K is l ess than the shift for the T at the same value of (3 since the J/ty is 
larger. The IP state does not shift, since it is not sensitive to perturbations at 
the origin. The change in the IP-IS splitting would then cause the derived a" 1 
to change upwards to 1.25(4) GeV if gluonic 0(a 2 ) effects were corrected. This 
is still within la of the original value. The expected shift in the ip' state is 0.005 
so the change in the 2S-1S splitting would be completely negligible compared to 
its statistical error. 

The value for a -1 is clearly different from that from T[16| or light hadron 



II spectroscopy at the same value of (3. In the quenched approximation we 
would expect > > a" 1 , reflecting the ordering of the momentum scales 
appropriate to the different quantities. In current results, the first inequality 



holds but the second one does not RIB]; this may reflect O(a) errors in present light 
hadron spectroscopy. Further calculations at different values of a will resolve this 
problem. 

The l D 2 state whose mass we have calculated is rather higher than that found 
for the ^(3770), thought to be a 3 Di state. From the spin splittings alone you 
would expect this difference. The ■0(3770) is also above threshold for decay to DD 
so quenching might have a significant effect on masses in this region, although the 
ratio of the width of the 0(3770) to its mass is still less than 1%. The 3 D\ has the 
same J PC quantum numbers as the 3 S\ and will appear as a third excited state 
in that channel. In order to observe such a state the cross-correlation between 
the meson correlators 3 S\ and the 3 D\ would have to be calculated and we have 
not attempted to do this here. 

Values for the wave function at the origin can be obtained as discussed in ref. 
@. If we include the (loc,loc) correlation function in a multi-exponential row 
fit we obtain a value of a 3 / 2 0(O) for the J/ty of 0.1535. This method does not 
yield a stable value for the excited states since the (loc, loc) correlation function 
does not distinguish different states very readily. A better method is take a 
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ratio of amplitudes from row and matrix fits ||. We use b(n sc ,m)/a(n sc ,m) and 
concentrate on the diagonal entries i.e. n sc = m =1 for J/ty and n sc = 2 for if}'. 
This gives a 3 / 2 ip(0) = 0.148(2) for J/tf and 0.13(1) for if;'. 

The leptonic width can be calculated from ^(0) using the Van-Royen Weis- 
skopf formula[|I7| at leading order. We obtain 5.4(5) keV for the J/^> in good 
agreement with the experimental value of 5.3(3) keV. The error we quote is dom- 
inated by the error in a -1 since this appears cubed. In principle we expect large 
corrections (~ 30%) to our value when a current correctly matched to the contin- 
uum current is included, instead of the naive lowest order current that we have 



used. We should apply small- components corrections to the current |L8| as well as 
a lattice-to-continuum renormalisation. The agreement with experiment should 
thus not be taken to be very significant at this stage. For the if)' the agreement 
with experiment is not so good. The experimental value is less than half of that 
for the IS and yet we obtain a ratio of 0.7 to the IS. This trend for excited states 
to have too large a value for ip(0) is again similar to that found in the T case. 
On improving the systematic error in our currents we would hope to notice an 
improvement here unless it is a feature of the quenched approximation. 

Spin splittings have been calculated for the ground S and P states. These are 
shown in detail in Figure (|2|). The agreement with experiment is good within 
expected systematic errors of 30 — 40 MeV. This would not be possible with- 
out tadpole-improvement of the spin-dependent terms. It was clear from the 
T spectrum || that splittings without tadpole improvement were about half the 
size of those with tadpole improvement. This would be an even bigger effect 
here where f3 and Uq are smaller. There is nevertheless some disagreement with 
experiment in Figure (^), and it is useful to find the source of this. There are 
sufficient experimental results for charmonium that the system provides a good 
test of the systematic removal of sources of error. 

From Table [11] we can see that the hyperfine splitting M( 3 S\) — M^Sq) has 
a very small statistical error. The difference from experiment then shows up 
clearly and is presumably a result of our systematic errors. There is again a 
30 — 40 MeV systematic error from higher order relativistic, discretisation and 
radiative corrections to the heavy quark action. This would be quite sufficient 
to explain the difference. Relativistic corrections are documented in JT|. The 
radiative corrections are 0(g 2 ) corrections to the coefficient of the a ■ B term 
beyond tadpole improvement. The discretisation errors are 0(a 2 ) errors in the B 
field and the hyperfine splitting is rather sensitive to these, as discussed below. 
We also expect quenching to have a significant effect, however. A comparison of 
T results on quenched and unquenched configurations shows an increase in the 
hyperfine splitting when unquenched (to 3 flavours) of between 30% and 50% 
|| [16[]. This can be explained largely on the basis of the difference between 
quenched and unquenched coupling constants and wavefunctions appearing in 
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the perturbative formula for the hyperfine splitting, 



AM hft =^p {mf . (16) 

For the J/ty case we might expect a similar shift of the hyperfine splitting on 
going to the full theory and this again would be sufficient to explain fully the 
deviation from experiment. One problem here is that the perturbative formulae 
are not quite as reliable as in the T case B. 



A calculation of the cc hyperfine splitting by the Fermilab group |19| gives a 
somewhat smaller value than ours. They use an improved Wilson fermion action 
for the heavy quarks and this approach has different systematic errors than ours. 

The case of the P state fine splittings is much more complicated, with an 
expected interplay of short and long range effects. In a potential model approach 
[ f2~0H two terms contribute - one proportional to (L ■ S) and the other proportional 
to (S12) where Si 2 = 4[3(s~i • n)(sa • n — s[ ■ s^)]. Here s\ and s 2 are the spins of 
the heavy quarks and n is an arbitrary unit vector. Evaluating these expectation 
values for 3 P Q , 3 Pi and 3 P 2 states of equal mass quarks allows us to compare ratios 
of the splittings, since the the expectation values of potentials that accompany 
these terms are the same for all P states. A useful ratio [ECU is 



r = M( X 2) - Mjxi) (u) 
M( Xl )-M( Xo ) { ] 

Experimental values are 0.48(1) for cc, 0.66(2) for bb (IP) and 0.58(3) for bb (2P). 
From a comparison of possible potentials to experiment the conventional picture 
emerges in which the spin-orbit potential appearing with L ■ S has both short 
and long range pieces, whereas the tensor potential appearing with Si 2 has only 
a short range piece. The short-range pieces can be related to 1-gluon exchange 
in perturbation theory and behave like 1/R 3 . The long-range piece comes from 
the scalar confining potential. Spin-dependent potentials can be extracted on 
the lattice from expectation values of Wilson loops with E and B field insertions 
along the time lines on either side. There it becomes clear that the 'same-side' 
spin-orbit potential is long-range, whereas the 'opposite-side' is short-range, as 
is the tensor potential |21 . 



We can extract values for the above ratio r of P spin splittings from our 
simulation and we find 1.2(2) for cc, clearly too large. The bb (IP) result at (3 — 
6.0 || is 0.7(3), which is consistent with experiment, but at (3 — 5.7 we obtain 
1.4(4) ||. It seems likely then that the disagreement with experiment arises 
from discretisation errors. At low j3 the predominant spin- dependent potential 
is the long-range spin-orbit piece, the shorter range pieces are not well resolved 
(compare [E2[ and [EjJ, for example). A pure L ■ S potential would give a value 



for r of 2 [p_0j (the pure tensor would give —0.4). In potential model language the 



long-range L ■ S term has undue dominance in our simulation. We also find that 
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the overall size of the P spin splittings, set by M(% 2 ) — M(xo), is too small. For 
bb at (3 = 6.0 this splitting was on the low side but in agreement with experiment 
within the error ||. For bb at (3 = 5.7 we obtain a result which is much too small 
0. Future calculations will concentrate on correcting discretisation errors to see 
if the results for charmonium at low (3 improve. 

Another possible discretisation error shows up in the fact that the centre of 
mass of the 3 P states comes out above the 1 P±. This happens both for this 
calculation and that of the T spectrum 0, but in both cases at a level within the 
expected systematic errors. One might expect, for example, that the hyperfine 
S ■ S interaction would contribute such a term to P states even in the absence of 
a wavefunction at the origin (see equation (16)) if the B field was smeared out 
over a plaquette as it is here. Experimental evidence so far indicates that there 



is no such splitting ||14j| , although it awaits confirmation. 

It seems likely that errors from the quenched approximation (and from dis- 
cretisation) are not so large for the P fine structure as for the S hyperfine splitting 
because the latter is determined by very short range phenomena. The hyperfine 
splitting can be thought of as resulting from delta-function S\ ■ S2 potential at 
the origin (see equation (16)), where S states have significant wave function. The 
quenched approximation causes larger effects at short distance scales because it 
appears, perturbatively, as an incorrect running of the coupling constant a(R) 
down to the origin from some R which is the important separation for quark and 
antiquark in the IP- IS splitting which is used to set the scale. P states have no 
wavefunction at the origin and in addition the short-range pieces of the relevant 
spin-dependent potentials have longer range than the delta function hyperfine 
potential. This should mean that the P fine structure can be determined accu- 
rately in a quenched calculation by a systematic improvement on this calculation, 
without having to unquench. 



5 Conclusions 

This represents a first calculation of the cc spectrum using NRQCD with spin- 
dependent terms. We include the leading relativistic and discretisation errors 
with tadpole- improved coefficients. We find a value of the lattice spacing from the 
IP-IS splitting which is different from that of the T on the same configurations 0. 
This is a clear indication of an effect from the quenched approximation. Another 
effect seen in the T spectrum itself, the difference in a -1 from the 2S-1S and 
IP-IS splittings, is not visible here above the statistical noise in the 2S state. 

With tadpole-improvement, the spin splittings agree with experiment at the 
level of the systematic error that we expect. The trend of these systematic 
errors is the same as that for the T spectrum and we would expect that, on 
including higher order terms, we could obtain better agreement. It seems likely 
that the major errors at present are discretisation effects and future calculations 
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will correct for these. One very good feature of the cc spectrum is that all the 
radial ground state S and P masses are known experimentally and so they can be 
used to gauge the effect of systematic improvement. Further calculations of the cc 
spectrum on lattices of different lattice spacing and on unquenched configurations 
would also provide useful checks of the systematic errors. A value for a s could 
be extracted from the IP-IS splitting in the same way that it was done using the 



T calculation || and a comparison with results from Wilson fermions p4| |25 
made. 

Calculations of the B c spectrum combining b and c propagators on these 
configurations will be reported shortly || . 
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